PRÁCTICA 1: APRENDIZAJE NO SUPERVISADO (DengAI) ¶
Nombre y apellidos: Cristian Guerrero Balber
Usuario VIU: cguerrerob@student.universidadviu.com
Esta actividad está dividida en dos grandes partes principales.
Ingeniería de características¶
Durante el desarrollo la primera parte, se va a analizar el dataset de DengAI de Driven data (dataset). Este análisis tendrá como objetivo entender y preprocesar mediante ingeniería de características el dataset con el fin de prepararlo para realizar una buena clusterización de datos, que será el objetivo final de la práctica.
La ingeniería de característica contendrá el estudio y tratamiento de errores en los datos, como el tratamiento de nulos o outliers.
Además, se hara uso de herramientas estadísticas y herramientas de asociación de características como la matriz de correlaciones, la matriz de similitud o algoritmos jerárquicos.
Clusterización¶
En esta segunda parte de la práctica se explorarán las distintas soluciones usando 4 algoritmos:
- K-MEANS
- DBSCAN
- GAUSSIAN MIXTURE MODELS (EXPECTATION-MAXIMIZATION)
- AFFINITY PROPAGATION
Se evaluarán los resultados, en su mayoría, utilizando la métrica de la silueta, ya que nuestros datos no están etiquetados.
Además, se utilizarán métricas exclusivas de cada modelo que no van a poder ser comparabes entre modelos pero ayudará a tomar una decisión final sobre qué modelo escoger y cuántos clústeres son los adecuados.
Finalmente, se llegará a una conclusión de porqué se selecciona el modelo que de un mejor rendimiento en base a las pruebas realizadas.
# Imports generales
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
import scipy.stats as stats
import io
from google.colab import files
seed = 42 # Semilla aleatoria arbitraria y constante a incluir en los algoritmos estocásticos para que los experimentos sean siempre reproducibles por el profesor.
# OJO: En los experimentos estocásticos que requieran varias iteraciones con distintas semillas, podéis incorporarla como seed+1, seed+2, etc.
def upload_files (index_fields=None):
uploaded = files.upload()
for fn in uploaded.keys():
print('User uploaded file "{name}" with length {length} bytes'.format(
name=fn, length=len(uploaded[fn])))
df = pd.read_csv(io.StringIO(uploaded[fn].decode('utf-8')), index_col = index_fields)
return df
# Subir el conjunto de entrenamiento sin variable objetivo (dengue_features_train.csv)
train = upload_files()
print(train.shape)
Saving dengue_features_train.csv to dengue_features_train.csv User uploaded file "dengue_features_train.csv" with length 287139 bytes (1456, 24)
train.head()
| city | year | weekofyear | week_start_date | ndvi_ne | ndvi_nw | ndvi_se | ndvi_sw | precipitation_amt_mm | reanalysis_air_temp_k | ... | reanalysis_precip_amt_kg_per_m2 | reanalysis_relative_humidity_percent | reanalysis_sat_precip_amt_mm | reanalysis_specific_humidity_g_per_kg | reanalysis_tdtr_k | station_avg_temp_c | station_diur_temp_rng_c | station_max_temp_c | station_min_temp_c | station_precip_mm | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | sj | 1990 | 18 | 1990-04-30 | 0.122600 | 0.103725 | 0.198483 | 0.177617 | 12.42 | 297.572857 | ... | 32.00 | 73.365714 | 12.42 | 14.012857 | 2.628571 | 25.442857 | 6.900000 | 29.4 | 20.0 | 16.0 |
| 1 | sj | 1990 | 19 | 1990-05-07 | 0.169900 | 0.142175 | 0.162357 | 0.155486 | 22.82 | 298.211429 | ... | 17.94 | 77.368571 | 22.82 | 15.372857 | 2.371429 | 26.714286 | 6.371429 | 31.7 | 22.2 | 8.6 |
| 2 | sj | 1990 | 20 | 1990-05-14 | 0.032250 | 0.172967 | 0.157200 | 0.170843 | 34.54 | 298.781429 | ... | 26.10 | 82.052857 | 34.54 | 16.848571 | 2.300000 | 26.714286 | 6.485714 | 32.2 | 22.8 | 41.4 |
| 3 | sj | 1990 | 21 | 1990-05-21 | 0.128633 | 0.245067 | 0.227557 | 0.235886 | 15.36 | 298.987143 | ... | 13.90 | 80.337143 | 15.36 | 16.672857 | 2.428571 | 27.471429 | 6.771429 | 33.3 | 23.3 | 4.0 |
| 4 | sj | 1990 | 22 | 1990-05-28 | 0.196200 | 0.262200 | 0.251200 | 0.247340 | 7.52 | 299.518571 | ... | 12.20 | 80.460000 | 7.52 | 17.210000 | 3.014286 | 28.942857 | 9.371429 | 35.0 | 23.9 | 5.8 |
5 rows × 24 columns
train.tail()
| city | year | weekofyear | week_start_date | ndvi_ne | ndvi_nw | ndvi_se | ndvi_sw | precipitation_amt_mm | reanalysis_air_temp_k | ... | reanalysis_precip_amt_kg_per_m2 | reanalysis_relative_humidity_percent | reanalysis_sat_precip_amt_mm | reanalysis_specific_humidity_g_per_kg | reanalysis_tdtr_k | station_avg_temp_c | station_diur_temp_rng_c | station_max_temp_c | station_min_temp_c | station_precip_mm | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1451 | iq | 2010 | 21 | 2010-05-28 | 0.342750 | 0.318900 | 0.256343 | 0.292514 | 55.30 | 299.334286 | ... | 45.00 | 88.765714 | 55.30 | 18.485714 | 9.800000 | 28.633333 | 11.933333 | 35.4 | 22.4 | 27.0 |
| 1452 | iq | 2010 | 22 | 2010-06-04 | 0.160157 | 0.160371 | 0.136043 | 0.225657 | 86.47 | 298.330000 | ... | 207.10 | 91.600000 | 86.47 | 18.070000 | 7.471429 | 27.433333 | 10.500000 | 34.7 | 21.7 | 36.6 |
| 1453 | iq | 2010 | 23 | 2010-06-11 | 0.247057 | 0.146057 | 0.250357 | 0.233714 | 58.94 | 296.598571 | ... | 50.60 | 94.280000 | 58.94 | 17.008571 | 7.500000 | 24.400000 | 6.900000 | 32.2 | 19.2 | 7.4 |
| 1454 | iq | 2010 | 24 | 2010-06-18 | 0.333914 | 0.245771 | 0.278886 | 0.325486 | 59.67 | 296.345714 | ... | 62.33 | 94.660000 | 59.67 | 16.815714 | 7.871429 | 25.433333 | 8.733333 | 31.2 | 21.0 | 16.0 |
| 1455 | iq | 2010 | 25 | 2010-06-25 | 0.298186 | 0.232971 | 0.274214 | 0.315757 | 63.22 | 298.097143 | ... | 36.90 | 89.082857 | 63.22 | 17.355714 | 11.014286 | 27.475000 | 9.900000 | 33.7 | 22.2 | 20.4 |
5 rows × 24 columns
Tras una primera visualización del dataset y como primer acercamiento al problema de clustering, se puede observar que existen variables temporales:
- year
- weekofyear
- week_start_date
Sería interesante ver la evolución temporal de todas las características a lo largo de los años. Se puede observar en la tabla que las fechas están ordenadas, por lo que simplemente se realizará una agregación y se obtendrá la media de los resultados.
plt.style.use('seaborn-v0_8-pastel')
fig, axs = plt.subplots(10, 2, figsize = (10, 15))
i, j = (0, 0)
for col in train.iloc[:, 4:].columns:
agg = train.groupby('year')[col].mean()
axs[i, j].plot(agg)
axs[i, j].set_title(col)
axs[i, j].tick_params(axis='x', labelsize=8)
if j < 1:
j += 1
else:
j = 0
i += 1
fig.suptitle('Evolución media de las características entre 1990 y 2010')
plt.tight_layout()
plt.show()
Se observa como la mayoria de las variables tienen una tendencia alcista a lo largo de los años, como reanalysis_max_air_temp_k. Otras, como reanalysis_min_air_temp_k, tienen una tendencia decreciente. Observando este comportamiento, se llega a la conclusión de que cada vez se alcanzan temperaturas más extremas según este dataset. Esto se confirma observando la variable reanalysis_tdtr_k, la cual tiene un salto creciente brusco a partir de 1999 aproximadamente. Esto es el rango máximo de temperatura en un día.
Observemos ahora la evolución temporal media en un año usando la columna weekofyear.
fig, axs = plt.subplots(10, 2, figsize = (10, 15))
i, j = (0, 0)
for col in train.iloc[:, 4:].columns:
agg = train.groupby('weekofyear')[col].mean()
axs[i, j].plot(agg)
axs[i, j].set_title(col)
axs[i, j].tick_params(axis='x', labelsize=8)
if j < 1:
j += 1
else:
j = 0
i += 1
fig.suptitle('Evolución media de las características anual')
plt.tight_layout()
plt.show()
Claramente la mayoría de las gráficas muestran un comportamiento estacional, lo cual tiene sentido debido a que las condiciones climáticas tienen mucho que ver con el momento del año en que se miden.
Por ejemplo, las máximas temperaturas se alcanzan sobre la semana 35, que coincide con estaciones de verano en la zona norte del planeta, mientras que las temperaturas más bajas son o al final del año o al principio, que corresponde al invierno.
Un detalle que me llama la atención es que, para algunas características como reanalysis_min_air_temp_k o reanalysis_relative_humidity_percent, existe un salto entre la semana 53 y la semana 1. Este salto no debería darse debido a que son períodos "continuos" de tiempo. Realmente, no es continuo ya que se ha realizado una agrupación y una media de los valores, pero al menos, debería de ser aproximado.
La pregunta que toca hacerse es, ¿son estas variables adecuadas para un trabajo de agrupación?
La variables year y week_start_date muestran una evolución temporal abierta, es decir, se podría extender hasta el infinito si existiesen datos, por lo tanto no es conveniente utilizarlas para una agrupación. La mayoría de modelos de clustering intentarían agrupar por cercanía de valores (1990 seria un grupo, y cerca de 1991, y asi sucesivamente). Esto no supone una ventaja para ver posibles agrupaciones útiles.
Sin embargo, weekofyear al ser una variable repetitiva en el tiempo, sí puede ser una buena variable agrupadora, ya que hay ciertas características subyacentes según en qué momento del año nos encontremos.
No obstante, sería conveniente realizar para esta última variable un estudio de la cantidad de instancias recogidas para cada valor, es decir, para cada semana del año.
train['weekofyear'].value_counts().plot(kind = 'bar')
<Axes: xlabel='weekofyear'>
Se puede observar que existen solamente 5 instancias de la semana 53. Esto puede explicar la falta de conexión entre la primera y última semana del año en algunas de las variables: la media de los datos de la semana 53 no es representativa.
Por todo lo anterior, la estrategia será la siguiente:
- Se eliminarán las características year y week_start_date.
- Se mantendrá la característica weekofyear, pero se eliminarán las instancias de la semana 53.
train.drop(["year", "week_start_date"], axis = 1, inplace = True)
train = train[train['weekofyear'] < 53]
Aclaración del significado de las variables¶
Tras esta primera limpieza basada en el sentido común y en los objetivos del problema, como ingeniero de datos he de realizar una buena comprensión de las características previamente al estudio estadístico para saber de qué estoy hablando.
Tras haber leído en Problem Description la definición de cada una, esto no es suficiente para comprender qué son algunas de ellas.
Por lo tanto, añado una breve literatura sobre algunas de ellas que son poco autodescriptivas.
Esto ayudará a entender mejor el dataset y, por lo tanto, tomar mejores decisiones en la construcción del modelo.
Características NDVI: Índice De Vegetación De Diferencia Normalizada¶
Explicado de la forma más sencilla posible, el Índice de Vegetación de Diferencia Normalizada (NDVI) mide el verdor y la densidad de la vegetación captada en una imagen de satélite. La vegetación sana tiene una curva de reflectancia espectral muy característica de la que podemos sacar partido calculando la diferencia entre dos bandas: la del rojo visible y la del infrarrojo cercano. El NDVI es esa diferencia expresada numéricamente entre -1 y 1 .
[...]
La caída en los valores también puede corresponder a cambios normales, como el momento de la cosecha, por lo que el NDVI debe contrastarse con otros datos disponibles.
[...]
Este índice está definido por valores que van de -1.0 a 1.0, donde los valores negativos están formados principalmente por nubes, agua y nieve, y los valores negativos cercanos a cero están formados principalmente por rocas y suelo descubierto.
Los valores muy pequeños (0,1 o menos) de la función NDVI corresponden a áreas sin rocas, arena o nieve.
Los valores moderados (de 0,2 a 0,3) representan arbustos y praderas, mientras que los valores grandes (de 0,6 a 0,8) indican bosques templados y tropicales.
EOSDA Crop Monitoring utiliza con éxito esta escala para mostrar a los agricultores qué partes de sus campos tienen una vegetación densa, moderada o escasa en un momento dado.
Además, según chatgpt, una resolución espacial de 0.5x0.5 significa que la superficie terrestre está dividida en cuadrículas de ese tamaño. En el ecuador, cada celda de esta cuadrícula corresponde aproximadamente a un área de 55 km x 55 km. Cada celda de la cuadrícula contendrá un valor promedio de NDVI basado en los datos satelitales para esa región.
Ahora se extraen las características mencionadas y se llama a la función.
NOAA's NCEP Climate Forecast System Reanalysis measurements, NOAA's GHCN daily climate data weather station measurements y¶
Climate Forecast System Reanalysis
El CFSR es un producto de reanálisis de tercera generación. Es un sistema global de alta resolución, acoplado entre la atmósfera, el océano, la superficie terrestre y el hielo marino, diseñado para proporcionar la mejor estimación del estado de estos dominios acoplados durante este periodo (NCAR).
Daily climate data weather station measurements Como el principal conjunto de datos de estaciones del NOAA, el Global Historical Climatology Network (GHCN) es una base de datos integrada de resúmenes climáticos diarios de estaciones de superficie terrestre de todo el mundo, e incluye un conjunto común de revisiones de control de calidad. Los Centros Nacionales de Información Ambiental (NCEI) del NOAA producen el conjunto de datos GHCN.
El algoritmo PERSIANN es un modelo de red neuronal artificial (ANN) basado en una red neuronal de retroalimentación multicapa (MFN) conocida como la Red de Propagación Contraria Modificada (Hsu, 1996). Este modelo híbrido consta de dos procesos. Primero, las imágenes infrarrojas (10.2–11.2 µm) se transforman en la capa oculta a través de un proceso de agrupamiento automático para formar lo que se conoce como un mapa de características autoorganizado (SOFM). El propósito de este proceso es detectar y clasificar patrones en los datos de entrada. (hess.copernicus)
Es decir, todas las variables que comienzan con reanalysis y precipitation_atm_mm (PERSIANN)son estimaciones indirectas por históricos y usando un modelo, mientras que las que comienzan con station son mediciones directas y, por lo tanto, más fiables.
Esta diferenciación es crucial si se quiere reducir dimensionalidad, ya que en caso de encontrar una alta correlación lineal entre dos variables: se optará por mantener las "station".
Una posible candidata a desaparecer a simple vista podría ser reanalysis_tdtr_k, que mide el rango de temperatura diurno, y existe otra variable llamada station_diur_temp_rng_c que mide exactamente lo mismo pero de forma directa.
Codificación de características categóricas¶
La variable city contiene la ciudad de donde se han tomado/calculado los datos. Sin embargo, para poder realizar cálculos con ella, se debe codificar a un valor númerico.
Como solo existen dos ciudades en el dataset, se codificarán como 0 y 1.
from sklearn.preprocessing import LabelEncoder
encoder = LabelEncoder()
train['city'] = encoder.fit_transform(train['city'])
<ipython-input-12-07c71286bdb0>:4: SettingWithCopyWarning: A value is trying to be set on a copy of a slice from a DataFrame. Try using .loc[row_indexer,col_indexer] = value instead See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy train['city'] = encoder.fit_transform(train['city'])
Correlación de características¶
Como apoyo visual de características es interesante realizar un análisis bivariable. Se expondrá un pairplot de todo el conjunto (análisis bivariante completo). Además, de esta forma se puede visualizar la forma de la distribución de cada una de ellas por separado (análisis univariante).
sns.pairplot(train, diag_kind = 'kde')
plt.show()
Debido a la gran cantidad de características que componen este dataset, es dificil discernir las relaciones entre cada una de ellas y saber qué característica tiene relación con cual.
Sin embargo, sí que se pueden observar rápidamente ciertos pares que tienen un aspecto bastante lineal, mientras la mayoría son más dispersas.
Como primer acercamiento a este análisis bivariante, se va a desglosar este pairplot en grupos de características asociadas según de acuerdo con el conocimiento del dominio que se ha adquirido en la descripción de variables. Es decir, se va a agrupar por tipos.
- Variables ndvi
- Variables reanalysis
- Variables station
train.columns
Index(['city', 'weekofyear', 'ndvi_ne', 'ndvi_nw', 'ndvi_se', 'ndvi_sw',
'precipitation_amt_mm', 'reanalysis_air_temp_k',
'reanalysis_avg_temp_k', 'reanalysis_dew_point_temp_k',
'reanalysis_max_air_temp_k', 'reanalysis_min_air_temp_k',
'reanalysis_precip_amt_kg_per_m2',
'reanalysis_relative_humidity_percent', 'reanalysis_sat_precip_amt_mm',
'reanalysis_specific_humidity_g_per_kg', 'reanalysis_tdtr_k',
'station_avg_temp_c', 'station_diur_temp_rng_c', 'station_max_temp_c',
'station_min_temp_c', 'station_precip_mm'],
dtype='object')
sns.pairplot(train[['ndvi_ne', 'ndvi_nw', 'ndvi_se', 'ndvi_sw']], diag_kind = 'kde')
plt.show()
Las variables ndvi están claramente correlacionadas, pero con cierta variabilidad también. Esto hace que sean unas variables interesantes seguramente a mantener.
sns.pairplot(train[['reanalysis_avg_temp_k', 'reanalysis_dew_point_temp_k',
'reanalysis_max_air_temp_k', 'reanalysis_min_air_temp_k',
'reanalysis_precip_amt_kg_per_m2',
'reanalysis_relative_humidity_percent', 'reanalysis_sat_precip_amt_mm',
'reanalysis_specific_humidity_g_per_kg', 'reanalysis_tdtr_k']], diag_kind = 'kde')
plt.show()
Claramente se visualizan dos variables altamente correlacionadas:
- reanalysis_humidity_g_per_kg con reanalysis_dew_point_temp_k
Además, se observan agrupaciones claras entre los datos de ciertos pares de variables, como reanalysis_relative_humidity_percent y renalysis_tdtr_k.
En cuanto a las distribuciones, en esta categoría se observan variables con un gran sesgo o incluso con más de una gaussiana en su distribución.
sns.pairplot(train[['station_avg_temp_c', 'station_diur_temp_rng_c', 'station_max_temp_c',
'station_min_temp_c', 'station_precip_mm']], diag_kind = 'kde')
plt.show()
Por último, las variables de tipo station, muestran cierta correlación entre algunas de ellas, pero no demasiada.
Ahora que se entiende mucho mejor qué es lo que representan las características mencionadas y como son las relaciones que existen entre los subgrupos de estas, se realizará un análisis de correlaciones para evitar multicolinealidad.
def plot_heatmap(df, labels, max = 1, min = -1):
#Se genera una máscara por el triángulo superior
mask = np.triu(np.ones_like(df, dtype=bool))
#Se establece la figura de matplotlib
f, ax = plt.subplots(figsize=(11, 11))
#Se genera un mapa de color divergente
cmap = sns.diverging_palette(230, 20, as_cmap=True)
#Se plotea el mapa de calor con la máscara y un ratio correcto
sns.heatmap(df, mask=mask, cmap=cmap, vmax=max, vmin = min, center=0,
xticklabels = labels, yticklabels = labels,
square=True, linewidths=.5, cbar_kws={"shrink": .5},
annot = True, fmt=".2f", annot_kws = {"size": 8})
plt.xticks(fontsize = 8)
plt.yticks(fontsize = 8)
#Código adaptado de https://seaborn.pydata.org/examples/many_pairwise_correlations.html
#Se genera un dataframe de correlaciones
train_corr = train.corr()
labels = train.columns.tolist()
plot_heatmap(train_corr, labels)
Observaciones de la matriz de correlaciones:
- Todo el grupo de características ndvi están claramente correlacionadas entre ellas, pero ninguna supera el valor de 0.9 en el coeficiente de Pearson y preferiblemente no se elimina ninguna.
- reanalysis_avg_temp_k y reanalysis_air_temp_k tienen una correlación de 0.9. Sin embargo, se mantendrán ambas por representar información lexicamente distinta (no es lo mismo una temperatura puntual que una media de la temperatura en ese día, aunque estén muy correlacionadas).
- Como era de esperar, reanalysis_tdtr_k tiene una alta relación con station_diur_temp_rng_c: su coeficiente de Pearson es de 0.88. Además, dicha variable tiene una correlación lineal de 0.92 con reanalysis_max_air_temp_k, por lo que su información puede quedar explicada con esta variable.En este caso y, aunque no supera el 0.9, se podría decidir eliminar la variable reanalysis_tdtr_k para evitar redundancias. Sin embargo, se tomará la decisión tras graficar próximamente estas tres variables y ver como quedan distribuidos sus datos.
- Aunque existen otras variables que presuntamente miden / calculan lo mismo como reanalysis_max_air_temp_k y station_max_temp_c, sus **c. de Pearson* no superan en ningún caso el 0.8. Por lo tanto, se mantienen.
- Existen ciertas variables que tienen una correlación lineal absoluta (c. de Pearson = 1). Sin lugar a dudas, al menos una de cada característica de las correlacionadas debe eliminarse del dataset. De otra manera, se produciría la multicolinealidad y causaría la inestabilidad de los coeficientes así como una reducción de la precisión.Estas son:
- reanalysis_sat_precip_amt_mm y precipitation_amt_mm. (se elimina esta última).
- reanalysis_specific_humidity_g_per_kg y reanalysis_dew_point_temp_k (se elimina esta última)
Se procede a mostrar las relaciones entre las características mencionadas.
sns.pairplot(train[['reanalysis_sat_precip_amt_mm', 'precipitation_amt_mm']],
diag_kind = 'kde')
plt.show()
Claramente, se eliminará una de ellas: precipitation_amt_mm
sns.pairplot(train[['reanalysis_specific_humidity_g_per_kg',
'reanalysis_dew_point_temp_k']],
diag_kind = 'kde')
plt.show()
Claramente, se eliminará una de ellas: reanalysis_dew_point_temp_k
sns.pairplot(train[['reanalysis_tdtr_k',
'station_diur_temp_rng_c',
'reanalysis_max_air_temp_k']],
diag_kind = 'kde')
plt.show()
Aquí no está tan clara la eliminación. Aunque es cierto que los c. de Pearson son altos y existe correlación, no es tan clara como los casos anteriores, existiendo una mayor dispersión.
Por lo tanto, se decide dejar estas variables.
train.drop(['precipitation_amt_mm', 'reanalysis_dew_point_temp_k'], axis = 1, inplace = True)
Estadísticas generales¶
Tras hacer una primera pasada de limpieza de características con una reducción de dimensionalidad, se procede a describir estadísticamente el dataset restante.
train.describe()
| city | weekofyear | ndvi_ne | ndvi_nw | ndvi_se | ndvi_sw | reanalysis_air_temp_k | reanalysis_avg_temp_k | reanalysis_max_air_temp_k | reanalysis_min_air_temp_k | reanalysis_precip_amt_kg_per_m2 | reanalysis_relative_humidity_percent | reanalysis_sat_precip_amt_mm | reanalysis_specific_humidity_g_per_kg | reanalysis_tdtr_k | station_avg_temp_c | station_diur_temp_rng_c | station_max_temp_c | station_min_temp_c | station_precip_mm | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| count | 1451.000000 | 1451.000000 | 1262.000000 | 1404.000000 | 1434.000000 | 1434.000000 | 1446.000000 | 1446.000000 | 1446.000000 | 1446.000000 | 1446.000000 | 1446.000000 | 1443.000000 | 1446.000000 | 1446.000000 | 1413.000000 | 1413.000000 | 1436.000000 | 1442.000000 | 1434.000000 |
| mean | 0.643005 | 26.412130 | 0.142294 | 0.130553 | 0.203783 | 0.202305 | 298.701852 | 299.225578 | 303.427109 | 295.719156 | 40.151819 | 82.161959 | 45.760388 | 16.746427 | 4.903754 | 27.185783 | 8.059328 | 32.452437 | 22.102150 | 39.326360 |
| std | 0.479279 | 14.964361 | 0.140531 | 0.119999 | 0.073860 | 0.083903 | 1.362420 | 1.261715 | 3.234601 | 2.565364 | 43.434399 | 7.153897 | 43.715537 | 1.542494 | 3.546445 | 1.292347 | 2.128568 | 1.959318 | 1.574066 | 47.455314 |
| min | 0.000000 | 1.000000 | -0.406250 | -0.456100 | -0.015533 | -0.063457 | 294.635714 | 294.892857 | 297.800000 | 286.900000 | 0.000000 | 57.787143 | 0.000000 | 11.715714 | 1.357143 | 21.400000 | 4.528571 | 26.700000 | 14.700000 | 0.000000 |
| 25% | 0.000000 | 13.000000 | 0.044950 | 0.049217 | 0.155087 | 0.144209 | 297.658929 | 298.257143 | 301.000000 | 293.900000 | 13.055000 | 77.177143 | 9.800000 | 15.557143 | 2.328571 | 26.300000 | 6.514286 | 31.100000 | 21.100000 | 8.700000 |
| 50% | 1.000000 | 26.000000 | 0.128817 | 0.121429 | 0.196050 | 0.189450 | 298.646429 | 299.289286 | 302.400000 | 296.200000 | 27.245000 | 80.301429 | 38.340000 | 17.087143 | 2.857143 | 27.414286 | 7.300000 | 32.800000 | 22.200000 | 23.850000 |
| 75% | 1.000000 | 39.000000 | 0.248483 | 0.216600 | 0.248846 | 0.246982 | 299.833571 | 300.207143 | 305.500000 | 297.900000 | 52.200000 | 86.357857 | 70.235000 | 17.978214 | 7.625000 | 28.157143 | 9.566667 | 33.900000 | 23.300000 | 53.900000 |
| max | 1.000000 | 52.000000 | 0.508357 | 0.454429 | 0.538314 | 0.546017 | 302.200000 | 302.928571 | 314.000000 | 299.900000 | 570.500000 | 98.610000 | 390.600000 | 20.461429 | 16.028571 | 30.800000 | 15.800000 | 42.200000 | 25.600000 | 543.300000 |
train.shape
(1451, 20)
Tratamiento de nulos¶
Gracias a la combinación de las funciones mostradas, se comprueba que no hay ninguna característica que no tenga valores nulos. Esto es un problema y hay que tratarlos.
Existen diferentes enfoques para tratar los valores nulos. Dependiendo de la naturaleza de las características y de la cantidad de valores faltantes que haya, se tomará la estrategia adecuada.
Por lo tanto, primero se va a calcular el % de valores nulos por característica:
def porc_nulls(train):
for col in train.columns:
nulos = len(train[train[col].isnull()])
porc_nulos = round((nulos / train.shape[0]) * 100, 0)
print(f"Valores nulos en {col}: {porc_nulos} %")
porc_nulls(train)
#La función ha sido creada porque se llamará de nuevo más adelante, así evito duplicar código
Valores nulos en city: 0.0 % Valores nulos en weekofyear: 0.0 % Valores nulos en ndvi_ne: 13.0 % Valores nulos en ndvi_nw: 3.0 % Valores nulos en ndvi_se: 1.0 % Valores nulos en ndvi_sw: 1.0 % Valores nulos en reanalysis_air_temp_k: 0.0 % Valores nulos en reanalysis_avg_temp_k: 0.0 % Valores nulos en reanalysis_max_air_temp_k: 0.0 % Valores nulos en reanalysis_min_air_temp_k: 0.0 % Valores nulos en reanalysis_precip_amt_kg_per_m2: 0.0 % Valores nulos en reanalysis_relative_humidity_percent: 0.0 % Valores nulos en reanalysis_sat_precip_amt_mm: 1.0 % Valores nulos en reanalysis_specific_humidity_g_per_kg: 0.0 % Valores nulos en reanalysis_tdtr_k: 0.0 % Valores nulos en station_avg_temp_c: 3.0 % Valores nulos en station_diur_temp_rng_c: 3.0 % Valores nulos en station_max_temp_c: 1.0 % Valores nulos en station_min_temp_c: 1.0 % Valores nulos en station_precip_mm: 1.0 %
La mayoría tiene un número muy reducido de valores nulos (< 4%), excepto ndvi_ne, que es de un 13%. Para porcentajes mucho más altos se recomienda no usar la característica, ya que la creación de datos postizos siempre distorsionan (en mayor o en menor medida).
Como el próximo paso sería realizar cambios en los valores nulos y esto distorsina el dataset, con el fin de tener un menor impacto se van a excluir, si hubiese, aquellas instancias donde haya una gran cantidad de características nulas. Estas instancias no aportan absolutamente nada (todas las características nulas) o poco (muchas características nulas) al dataset e introducirán ruido en cambio.
Dicho de otro modo, cuanto más cantidad de características nulas tenga una instancia, más se parecería el rellenado de nulos a la adición de filas completamente nuevas, por inteligente que fuera dicho rellenado.
Primero, se comprueba si existe alguna instancia de este tipo:
'''Se genera un df con los booleanos True o False de na, es decir, cada valor na
será sustituido por True, sino, por False'''
train_na = train.isna()
#Se genera una serie para usar de filtro con tan solo las filas con algun na
filter_na = train.isna().any(axis = 1)
'''Ahora, se calcula el ratio entre columnas con na y el total de columnas.
Como criterio, se establecerá que, en caso de que haya un 50% o más de
características nulas, se eliminará la instancia.'''
ind_high_na = [] #Se crea una lista para almacenar indices de ins. a elim.
for i, row in train_na[filter_na].iterrows():
ratio_na = row.sum() / row.size
if ratio_na >= 0.5:
ind_high_na.append(i)
#Ahora, se eliminan dichos indices del dataset
train.drop(index = ind_high_na, inplace = True)
Se vuelven a calcular los porcentajes de nulos
porc_nulls(train)
Valores nulos en city: 0.0 % Valores nulos en weekofyear: 0.0 % Valores nulos en ndvi_ne: 13.0 % Valores nulos en ndvi_nw: 3.0 % Valores nulos en ndvi_se: 1.0 % Valores nulos en ndvi_sw: 1.0 % Valores nulos en reanalysis_air_temp_k: 0.0 % Valores nulos en reanalysis_avg_temp_k: 0.0 % Valores nulos en reanalysis_max_air_temp_k: 0.0 % Valores nulos en reanalysis_min_air_temp_k: 0.0 % Valores nulos en reanalysis_precip_amt_kg_per_m2: 0.0 % Valores nulos en reanalysis_relative_humidity_percent: 0.0 % Valores nulos en reanalysis_sat_precip_amt_mm: 0.0 % Valores nulos en reanalysis_specific_humidity_g_per_kg: 0.0 % Valores nulos en reanalysis_tdtr_k: 0.0 % Valores nulos en station_avg_temp_c: 2.0 % Valores nulos en station_diur_temp_rng_c: 2.0 % Valores nulos en station_max_temp_c: 1.0 % Valores nulos en station_min_temp_c: 0.0 % Valores nulos en station_precip_mm: 1.0 %
train.shape
(1446, 20)
Al haber eliminado aquellas instancias con un alto porcentaje de caracteríscicas nulas, nos quedamos con un dataset donde muchas de las características con balores bajos aún se han reducido más y tan solo una característica tiene valores nulos: ndvi_ne con un 13%.
Además, no se han perdido demasiadas instancias, ya que hemos pasado de un dataset de 1456 filas a 1446 (solo se han perdido 10 instancias, que es asumible entre miles de ellas).
Al ser un 13% un valor relativamente alto de valores nulos (pero no tanto como para prescindir de la característica), se optará por una manera de cumplimentar los nulos un poco más sofisticada que la sustitución por la media o algun otro parámetro estadístico centralizador.
Se usará un modelo predictor basado en regresión lineal.
Sin embargo, para el resto de características con un porcentaje <3% de valores nulos, se cumplimentará con la media o la mediana dependiendo del sesgo de su distribución.
Para la regresión lineal, tiene sentido calcular el valor de ndvi_ne en función de las otras tres del mismo grupo (otras ndvi), ya que todas provienen del mismo sistema de medición y, como se comprobó anteriormente, todas tienen cierta correlación entre ellas.
En los modelos de regresión lineal, las características involucradas idealmente deberían tener una distribución normal, ya que este modelo es sensible a outliers.
Por lo tanto, primero se va a visualizar la forma de la distribución de dichas variables así como su comportamiento entre ellas a pares.
Análisis de la simetría¶
Para estudiar la simetría de estas variables más allá de lo visual, se va a usar el sesgo (skewness) como indicador estadístico (también conocido como tercer momento. Este sesgo mide la asimetría de la distribución.
Para evaluar el sesgo numéricamente y saber si la distribución es asimétrica o no, seguiremos las siguientes categorizaciones (datacamp):
- (-0.5, 0.5): bajo sesgo o aproximadamente simétrica.
- (-1, -0.5) U (0.5, 1): moderadamente sesgada
- (-inf, -1) U (1, inf): altamente sesgada
Además, se aprovechará este cálculo para el resto de características y su cumplimentación por la media o la mediana.
train.info()
<class 'pandas.core.frame.DataFrame'> Index: 1446 entries, 0 to 1455 Data columns (total 20 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 city 1446 non-null int64 1 weekofyear 1446 non-null int64 2 ndvi_ne 1257 non-null float64 3 ndvi_nw 1399 non-null float64 4 ndvi_se 1429 non-null float64 5 ndvi_sw 1429 non-null float64 6 reanalysis_air_temp_k 1446 non-null float64 7 reanalysis_avg_temp_k 1446 non-null float64 8 reanalysis_max_air_temp_k 1446 non-null float64 9 reanalysis_min_air_temp_k 1446 non-null float64 10 reanalysis_precip_amt_kg_per_m2 1446 non-null float64 11 reanalysis_relative_humidity_percent 1446 non-null float64 12 reanalysis_sat_precip_amt_mm 1443 non-null float64 13 reanalysis_specific_humidity_g_per_kg 1446 non-null float64 14 reanalysis_tdtr_k 1446 non-null float64 15 station_avg_temp_c 1413 non-null float64 16 station_diur_temp_rng_c 1413 non-null float64 17 station_max_temp_c 1436 non-null float64 18 station_min_temp_c 1442 non-null float64 19 station_precip_mm 1434 non-null float64 dtypes: float64(18), int64(2) memory usage: 237.2 KB
sesgos = {} #En este diccionario se almacenarán los sesgos por característica
tipos_sesgos = ('Bajo sesgo o aproximadamente simétrica',
'Moderadamente sesgada',
'Altamente sesgada')
for col in train.columns:
sesgo = stats.skew(train[col], nan_policy = 'omit')
if -0.5 < sesgo < 0.5:
tipo = tipos_sesgos[0]
elif -1 < sesgo <= -0.5 or 0.5 <= sesgo < 1:
tipo = tipos_sesgos[1]
elif sesgo <= -1 or sesgo >= 1:
tipo = tipos_sesgos[2]
sesgos[col] = tipo
print(f"Característica: {col}; Sesgo: {round(sesgo, 2)}, {tipo}")
#Referencia: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.skew.html
Característica: city; Sesgo: -0.6, Moderadamente sesgada Característica: weekofyear; Sesgo: 0.0, Bajo sesgo o aproximadamente simétrica Característica: ndvi_ne; Sesgo: -0.1, Bajo sesgo o aproximadamente simétrica Característica: ndvi_nw; Sesgo: -0.01, Bajo sesgo o aproximadamente simétrica Característica: ndvi_se; Sesgo: 0.57, Moderadamente sesgada Característica: ndvi_sw; Sesgo: 0.75, Moderadamente sesgada Característica: reanalysis_air_temp_k; Sesgo: -0.08, Bajo sesgo o aproximadamente simétrica Característica: reanalysis_avg_temp_k; Sesgo: -0.19, Bajo sesgo o aproximadamente simétrica Característica: reanalysis_max_air_temp_k; Sesgo: 0.85, Moderadamente sesgada Característica: reanalysis_min_air_temp_k; Sesgo: -0.67, Moderadamente sesgada Característica: reanalysis_precip_amt_kg_per_m2; Sesgo: 3.38, Altamente sesgada Característica: reanalysis_relative_humidity_percent; Sesgo: 0.57, Moderadamente sesgada Característica: reanalysis_sat_precip_amt_mm; Sesgo: 1.74, Altamente sesgada Característica: reanalysis_specific_humidity_g_per_kg; Sesgo: -0.54, Moderadamente sesgada Característica: reanalysis_tdtr_k; Sesgo: 1.07, Altamente sesgada Característica: station_avg_temp_c; Sesgo: -0.57, Moderadamente sesgada Característica: station_diur_temp_rng_c; Sesgo: 0.84, Moderadamente sesgada Característica: station_max_temp_c; Sesgo: -0.26, Bajo sesgo o aproximadamente simétrica Característica: station_min_temp_c; Sesgo: -0.31, Bajo sesgo o aproximadamente simétrica Característica: station_precip_mm; Sesgo: 2.98, Altamente sesgada
Cumplimentación de características con bajo % de nulos¶
La estrategia será la siguiente:
- Si el sesgo es Bajo sesgo o aproximadamente simétrica, se cumplimentarán por la media.
- Si el sesgo es Moderadamente sesgada o Altamente sesgada, se cumplimentará por la mediana.
El motivo de esta diferenciación es que, por norma, la media es la mejor medida centralizadora estadística. Sin embargo, es muy susceptible a los outliers y, por lo tanto, a la asimetria.
La mediana en cambio es más robusta ante distribuciones asimétricas con colas debido a outliers, por lo que se convierte en la mejor opción para este tipo de características.
for col in train.columns:
if col == 'ndvi_ne':
continue
if sesgos[col] == tipos_sesgos[0]:
#train[col].fillna(train[col].mean(), inplace = True) #Obsoleto a partir de pandas 3.0, me apareció mensaje de error
train.fillna({col: train[col].mean()}, inplace = True)
else:
#train[col].fillna(train[col].median(), inplace = True) #Obsoleto a partir de pandas 3.0, me apareció mensaje de error
train.fillna({col: train[col].median()}, inplace = True)
porc_nulls(train)
Valores nulos en city: 0.0 % Valores nulos en weekofyear: 0.0 % Valores nulos en ndvi_ne: 13.0 % Valores nulos en ndvi_nw: 0.0 % Valores nulos en ndvi_se: 0.0 % Valores nulos en ndvi_sw: 0.0 % Valores nulos en reanalysis_air_temp_k: 0.0 % Valores nulos en reanalysis_avg_temp_k: 0.0 % Valores nulos en reanalysis_max_air_temp_k: 0.0 % Valores nulos en reanalysis_min_air_temp_k: 0.0 % Valores nulos en reanalysis_precip_amt_kg_per_m2: 0.0 % Valores nulos en reanalysis_relative_humidity_percent: 0.0 % Valores nulos en reanalysis_sat_precip_amt_mm: 0.0 % Valores nulos en reanalysis_specific_humidity_g_per_kg: 0.0 % Valores nulos en reanalysis_tdtr_k: 0.0 % Valores nulos en station_avg_temp_c: 0.0 % Valores nulos en station_diur_temp_rng_c: 0.0 % Valores nulos en station_max_temp_c: 0.0 % Valores nulos en station_min_temp_c: 0.0 % Valores nulos en station_precip_mm: 0.0 %
Se puede comprobar que, ahora, no existen valores nulos excepto para la característica, ndvi_ne. Queda el dataset preparado para la regresión lineal.
Cumplimentación de la característica con alto % de nulos: Modelo de Regresión Lineal¶
Según las observaciones y los cálculos, entre las variables ndvi no hay ninguna altamente sesgada y, por lo tanto, concluyo que son variables aceptables para su uso en una regresión lineal.
from sklearn.linear_model import LinearRegression
Ahora, se va a dividir el df en dos:
- df_to_predict: parte del dataframe que contiene las instancias con valores nulos de ndvi_ne. Se usará para predecir la variable en cuestión.
- df_subtrain: parte del dataframe que no contiene ninguna instancia con valores nulos. Se usará para realizar el entrenamiento que predecirá la variable ndvi_ne.
df_to_predict = train[train.isna().any(axis = 1)][['ndvi_ne', 'ndvi_nw', 'ndvi_se', 'ndvi_sw']]
df_subtrain = train[~train.isna().any(axis = 1)][['ndvi_ne', 'ndvi_nw', 'ndvi_se', 'ndvi_sw']]
Del dataset df_subtrain (sub- porque no es el dataset principal de entrenamiento de este ejercicio en general), se va a dividir las columnas en x_train, y
from sklearn.model_selection import train_test_split
#https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.train_test_split.html
X = df_subtrain[['ndvi_nw', 'ndvi_se', 'ndvi_sw']] #Variables independientes o explicativas
y = df_subtrain[['ndvi_ne']] #Variable dependiente o explicada
Se separa el dataset en los datos de entrenamiento y validación (aunque por convención se les llame X_test e y_test).
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.2,
random_state = 42,
shuffle = True)
Se entrena el modelo linear con los datos de entrenamiento
linear_model = LinearRegression()
linear_model.fit(X_train, y_train)
#https://scikit-learn.org/dev/modules/generated/sklearn.linear_model.LinearRegression.html
LinearRegression()In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
LinearRegression()
Se comprueba mediante el indicador r2 si la regresión se ajusta bien a los datos de test.
linear_model.score(X_test, y_test)
0.727325236857673
Un r2 de 0.73 es bastante aceptable. No obstante, y aunque el desarrollo seguido tiene bastante lógica, en ML es esencial hacer ciertas pruebas para comprobar que el modelo escogido es lo más óptimo posible.
Por lo tanto, se realizará una prueba extra. Se creará el mismo modelo pero sin excluir a ninguna de las variables del dataset (sin estudiar, siquiera, sus distribuciones).
Se realiza el mismo proceso que le anterior pero sin excuir ninguna variable.
df_to_predict = train[train.isna().any(axis = 1)]
df_subtrain = train[~train.isna().any(axis = 1)]
X = df_subtrain.drop('ndvi_ne', axis = 1) #Variables independientes o explicativas
y = df_subtrain[['ndvi_ne']] #Variable dependiente o explicada
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.2,
random_state = 42,
shuffle = True)
linear_model = LinearRegression()
linear_model.fit(X_train, y_train)
linear_model.score(X_test, y_test)
0.7399935557180093
Se ha obtenido una mejor puntuación r2 a la anterior, por lo tanto, este será el modelo que se usará para la predicción de valores faltantes.
En ocasiones, aunque no todas las características estén muy correlacionadas con una variable a explicar, muchas pequeñas variaciones a tener en cuenta de cada una (pequeñas correlaciones) pueden dar mejores resultados que buenas correlaciones de pocas características. Aquí reside la potencia del ML.
Ahora, se realiza el proceso de separación de variables independientes de la dependiente sobre el dataset con valores faltantes.
X = df_to_predict.drop('ndvi_ne', axis = 1)
Se predice, en base al modelo creado, los valores faltantes de ndvi_ne
y_pred = linear_model.predict(X)
Se sustituyen los valores predichos en el dataset.
dict_index_values = {}
j = 0
for i, _ in df_to_predict['ndvi_ne'].items():
dict_index_values[i] = y_pred[j]
j += 1
train.fillna({'ndvi_ne': dict_index_values}, inplace = True)
train['ndvi_ne'] = train['ndvi_ne'].astype(float)
porc_nulls(train)
Valores nulos en city: 0.0 % Valores nulos en weekofyear: 0.0 % Valores nulos en ndvi_ne: 0.0 % Valores nulos en ndvi_nw: 0.0 % Valores nulos en ndvi_se: 0.0 % Valores nulos en ndvi_sw: 0.0 % Valores nulos en reanalysis_air_temp_k: 0.0 % Valores nulos en reanalysis_avg_temp_k: 0.0 % Valores nulos en reanalysis_max_air_temp_k: 0.0 % Valores nulos en reanalysis_min_air_temp_k: 0.0 % Valores nulos en reanalysis_precip_amt_kg_per_m2: 0.0 % Valores nulos en reanalysis_relative_humidity_percent: 0.0 % Valores nulos en reanalysis_sat_precip_amt_mm: 0.0 % Valores nulos en reanalysis_specific_humidity_g_per_kg: 0.0 % Valores nulos en reanalysis_tdtr_k: 0.0 % Valores nulos en station_avg_temp_c: 0.0 % Valores nulos en station_diur_temp_rng_c: 0.0 % Valores nulos en station_max_temp_c: 0.0 % Valores nulos en station_min_temp_c: 0.0 % Valores nulos en station_precip_mm: 0.0 %
Con esto finaliza el tratamiento de valores nulos.
Tratamiento de duplicados¶
Una columna cuyos valores sean siempre los mismos, no suele aportar nada al modelo. Este es un problema típico en características categóricas. Sin embargo, el dataset que se está tratando, tiene todas las características numéricas continuas. Por lo tanto, no tiene sentido realizar este estudio por característica.
Sin embargo, se puede realizar un checkeo rápido de que, al menos, no haya instancias duplicadas. Esto podría darse por un error de tratamiento del dataset previo a este estudio o por pura coincidencia (lo cual es improbable dado el número de dimensiones que se están tratando).
train.duplicated().any()
False
No existe ninguna instancia completamente duplicada. Por lo tanto, se continua sin cambiar el dataset.
Outliers¶
Como se ha visto en el apartado de Tratamiento de nulos, existen ciertas características con un sesgo alto. La mayoría de algoritmos de ML, son sensibles a estos sesgos, por lo que es conveniente visualizar los outliers y, en caso de que tenga sentido, tratarlos.
sesgos
{'city': 'Moderadamente sesgada',
'weekofyear': 'Bajo sesgo o aproximadamente simétrica',
'ndvi_ne': 'Bajo sesgo o aproximadamente simétrica',
'ndvi_nw': 'Bajo sesgo o aproximadamente simétrica',
'ndvi_se': 'Moderadamente sesgada',
'ndvi_sw': 'Moderadamente sesgada',
'reanalysis_air_temp_k': 'Bajo sesgo o aproximadamente simétrica',
'reanalysis_avg_temp_k': 'Bajo sesgo o aproximadamente simétrica',
'reanalysis_max_air_temp_k': 'Moderadamente sesgada',
'reanalysis_min_air_temp_k': 'Moderadamente sesgada',
'reanalysis_precip_amt_kg_per_m2': 'Altamente sesgada',
'reanalysis_relative_humidity_percent': 'Moderadamente sesgada',
'reanalysis_sat_precip_amt_mm': 'Altamente sesgada',
'reanalysis_specific_humidity_g_per_kg': 'Moderadamente sesgada',
'reanalysis_tdtr_k': 'Altamente sesgada',
'station_avg_temp_c': 'Moderadamente sesgada',
'station_diur_temp_rng_c': 'Moderadamente sesgada',
'station_max_temp_c': 'Bajo sesgo o aproximadamente simétrica',
'station_min_temp_c': 'Bajo sesgo o aproximadamente simétrica',
'station_precip_mm': 'Altamente sesgada'}
Se crea un generador para incrementar las veces que hagan falta dos coordenadas que servirán para graficar n veces boxplot en función del criterio mencionado.
def increment_index(k):
for i in range (k):
for j in range (k):
yield i, j
#Referencia: https://ellibrodepython.com/yield-python
import math
def plotear_boxplots(tipo):
global train
global sesgos
cols_to_plot = []
for col, value in sesgos.items():
if value == tipo:
#plt.boxplot(train[col])
cols_to_plot.append(col)
'''Como quiero crear una serie de subplots y que se pinten en un grid
cuadrado, calculo la raíz cuadrada k (entera redondeada hacia arriba) de la
cantidad de gráficas a pintar. El resultado será una matriz de subplots de kxk'''
k = math.ceil(math.sqrt(len(cols_to_plot)))
#Se crea el generador para el valor de k como límite del grid
incr_i_j = increment_index(k)
#Se crea el grid
fig, axs = plt.subplots(k, k, figsize = (10, 10))
plt.suptitle(f'Boxplots de las características de tipo: {tipo}')
for col in cols_to_plot:
i, j = next(incr_i_j)
axs[i, j].boxplot(train[col], patch_artist = True,
boxprops={'facecolor': 'skyblue'},
medianprops={'color': 'blue', 'linewidth': 1.5,
'linestyle': '--'},
flierprops = {'markerfacecolor': 'red',
'markersize': 4,
'linestyle': 'none'})
axs[i, j].set_title(f'{col}', fontsize = 8)
axs[i, j].set_xticks([])
plt.tight_layout(rect = [0, 0, 1, 0.96])
plt.show()
#Referencias:
#https://matplotlib.org/stable/gallery/subplots_axes_and_figures/subplots_demo.html
#https://matplotlib.org/3.1.1/api/_as_gen/matplotlib.pyplot.boxplot.html
#chatgpt
for tipo in tipos_sesgos:
plotear_boxplots(tipo)
Se puede observar que los boxplots del tipo "Altamente Sesgada" presentan ciertos valores extremadamente más alejados de los "outliers normales", o mejor dicho, los no tan extremos.
Aunque los outliers, por definición, son valores atípicos, al haber un número relativamente elevado de ellos no se van a eliminar debido a que podrían formar su propio cluster o contener información relevante de sus instancias completas.
Sin embargo, sí que se va a tratar los outliers más alejados. Estos son muy pocos y no son representativos para la formación de clusteres o, para un futuro ejercicio, predicciones. Por el contrario, pueden dificultar la tarea para la mayoría de algoritmos.
En lugar de eliminar las instancias completas correspondientes a estos puntos extremadamente atípicos, se va a sustituir su valor por la media, que seguirá estándo desplazada hacia arriba (debido al sesgo) pero no será tan extrema.
train['reanalysis_precip_amt_kg_per_m2'] = train['reanalysis_precip_amt_kg_per_m2'].apply(
lambda x: x if x < 320 else train['reanalysis_precip_amt_kg_per_m2'].mean())
train['station_precip_mm'] = train['station_precip_mm'].apply(
lambda x: x if x < 320 else train['station_precip_mm'].mean())
train['reanalysis_sat_precip_amt_mm'] = train['reanalysis_sat_precip_amt_mm'].apply(
lambda x: x if x < 260 else train['reanalysis_sat_precip_amt_mm'].mean())
plotear_boxplots(tipos_sesgos[2])
Con esto, se finaliza el tratamiento de outliers.
Creación de nuevas características¶
Para crear nuevas características a partir de las existentes, es útil centrarse en la naturaleza de los datos y las relaciones que podrían influir en la segmentación de clústeres.
1. Promedio y Rango de NDVI¶
Un promedio de NDVI puede representar la salud general de la vegetación en la zona, mientras que el rango puede indicar la diversidad de la vegetación.
- Nuevas características:
ndvi_avg: Promedio dendvi_ne,ndvi_nw,ndvi_se,ndvi_sw.ndvi_range: Rango de NDVI (máximo - mínimo).
2. Índices de Temperaturas¶
Estos índices pueden reflejar las condiciones climáticas extremas y su variabilidad, que pueden influir en el tipo de vegetación y el uso del suelo.
- Nuevas características:
temp_range_c: Rango de temperatura (máxima - mínima) usandostation_max_temp_cystation_min_temp_c.total_avg_air_temp: Promedio de los distintos promedios de temperaturas del aire utilizandoreanalysis_air_temp_kystation_avg_temp_c.temp_range_k: Diferencia entrereanalysis_max_air_temp_kyreanalysis_min_air_temp_ka lo largo de las semanas del año.avg_temp_range_c: media de los rangos de temperaturas obtenidas de fuentes diferentes.
3. Precipitación media¶
La precipitación total es un factor clave que afecta la agricultura y el crecimiento de la vegetación. Una medida compuesta puede ser útil para evaluar la disponibilidad de agua en diferentes áreas.
- Nuevas características:
avg_precip: media dereanalysis_precip_amt_kg_per_m2ystation_precip_mm. De forma aproximada, un kg de agua por metro cuadrado equivale a una altura de 1 mm de agua, ya que este volumen equivale a un litro (y la densidad relativa del agua es 1). Por lo tanto, son interoperables sin necesidad de ninguna transformación.
Para identificar claramente a todas estas nuevas caraceterísticas, llevarán el prefijo calc_(de calculadas o calculated, ya que todas las variables las pongo en inglés).
# NDVIs
## Media ndvi
train['calc_ndvi_avg'] = train[['ndvi_ne', 'ndvi_nw',
'ndvi_se', 'ndvi_sw']].mean(axis = 1)
## Rangos ndvi
train['calc_ndvi_range'] = train[['ndvi_ne', 'ndvi_nw',
'ndvi_se', 'ndvi_sw']].max(axis = 1)-train[['ndvi_ne',
'ndvi_nw', 'ndvi_se', 'ndvi_sw']].min(axis = 1)
# Temperaturas
## temp_range
train['calc_temp_range_c'] = train['station_max_temp_c'] - train['station_min_temp_c']
## total_avg_air_temp
train['calc_total_avg_air_temp'] = (train['reanalysis_air_temp_k'] - 273.15 + train['station_avg_temp_c']) / 2
## temp_range_k
train['calc_temp_range_k'] = train['reanalysis_max_air_temp_k'] - train['reanalysis_min_air_temp_k']
## avg_temp_range
train['calc_avg_temp_range_c'] = (train['calc_temp_range_k'] - 273.15 + train['calc_temp_range_c']) / 2
# Precipitaciones
## avg_precip
train['calc_avg_precip'] = train[['reanalysis_precip_amt_kg_per_m2', 'station_precip_mm']].mean(axis=1)
#Referencias:
#https://pandas.pydata.org/docs/reference/api/pandas.DataFrame.max.html
train.columns
Index(['city', 'weekofyear', 'ndvi_ne', 'ndvi_nw', 'ndvi_se', 'ndvi_sw',
'reanalysis_air_temp_k', 'reanalysis_avg_temp_k',
'reanalysis_max_air_temp_k', 'reanalysis_min_air_temp_k',
'reanalysis_precip_amt_kg_per_m2',
'reanalysis_relative_humidity_percent', 'reanalysis_sat_precip_amt_mm',
'reanalysis_specific_humidity_g_per_kg', 'reanalysis_tdtr_k',
'station_avg_temp_c', 'station_diur_temp_rng_c', 'station_max_temp_c',
'station_min_temp_c', 'station_precip_mm', 'calc_ndvi_avg',
'calc_ndvi_range', 'calc_temp_range_c', 'calc_total_avg_air_temp',
'calc_temp_range_k', 'calc_avg_temp_range_c', 'calc_avg_precip'],
dtype='object')
Al haber introducido nuevas variables, veamos de nuevo rápidamente su matriz de corrleaciones.
train_corr = train.corr()
plot_heatmap(train_corr, train.columns.tolist())
#Código adaptado de https://seaborn.pydata.org/examples/many_pairwise_correlations.html
Salta a la vista que se han creado algunas variables demasiado correlacionadas. Esto puede crear problemas de redundancia en el modelo. Las dos más preocupantes y que se deben eliminar del modelo por tener un coeficiente de Pearson igual o mayor a 0.96 son:
calc_temp_range_k(a eliminar) conreanalysis_tdtr_k(0.97)calc_avg_temp_range_c(a eliminar) conreanalysis_tdtr_k(0.96) y concalc_temp_range_kprecisamente.
Se ha comprobado que estas variables no aportan explicación extra al modelo, por lo que serán eliminadas para evitar multicolinealidad.
train.drop(['calc_temp_range_k', 'calc_avg_temp_range_c'],
axis = 1, inplace = True)
Relación Jerárquica de Características¶
Aunque ya se ha realizado una primera transformación del dataset en la exploración de datos (como el estudio de correlaciones y la eliminación de características), en este apartado se va a realizar una ingeniería de características un poco más avanzada, viendo posibles agrupaciones entre características similares.
Para ello, se va a usar uno de los algoritmos de clasificación no supervisada: la clasisficación Jerárquica.
Para ello, como primer paso se transpondrá el dataset, para convertir las características a instancias y, como columnas, simplemente tendremos los índices (0, 1, 2...).
train_transposed = train.T
train_transposed.head()
| 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | ... | 1446 | 1447 | 1448 | 1449 | 1450 | 1451 | 1452 | 1453 | 1454 | 1455 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| city | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.0000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | ... | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 |
| weekofyear | 18.000000 | 19.000000 | 20.000000 | 21.000000 | 22.0000 | 23.000000 | 24.000000 | 25.000000 | 26.000000 | 27.000000 | ... | 16.000000 | 17.000000 | 18.000000 | 19.000000 | 20.000000 | 21.000000 | 22.000000 | 23.000000 | 24.000000 | 25.000000 |
| ndvi_ne | 0.122600 | 0.169900 | 0.032250 | 0.128633 | 0.1962 | 0.164221 | 0.112900 | 0.072500 | 0.102450 | 0.093376 | ... | 0.231486 | 0.239743 | 0.260814 | 0.168686 | 0.263071 | 0.342750 | 0.160157 | 0.247057 | 0.333914 | 0.298186 |
| ndvi_nw | 0.103725 | 0.142175 | 0.172967 | 0.245067 | 0.2622 | 0.174850 | 0.092800 | 0.072500 | 0.146175 | 0.121550 | ... | 0.294686 | 0.259271 | 0.255786 | 0.158500 | 0.272500 | 0.318900 | 0.160371 | 0.146057 | 0.245771 | 0.232971 |
| ndvi_se | 0.198483 | 0.162357 | 0.157200 | 0.227557 | 0.2512 | 0.254314 | 0.205071 | 0.151471 | 0.125571 | 0.160683 | ... | 0.331657 | 0.307786 | 0.257771 | 0.133071 | 0.258271 | 0.256343 | 0.136043 | 0.250357 | 0.278886 | 0.274214 |
5 rows × 1446 columns
Ahora, como índices tendremos las características.
names = train_transposed.index
names
Index(['city', 'weekofyear', 'ndvi_ne', 'ndvi_nw', 'ndvi_se', 'ndvi_sw',
'reanalysis_air_temp_k', 'reanalysis_avg_temp_k',
'reanalysis_max_air_temp_k', 'reanalysis_min_air_temp_k',
'reanalysis_precip_amt_kg_per_m2',
'reanalysis_relative_humidity_percent', 'reanalysis_sat_precip_amt_mm',
'reanalysis_specific_humidity_g_per_kg', 'reanalysis_tdtr_k',
'station_avg_temp_c', 'station_diur_temp_rng_c', 'station_max_temp_c',
'station_min_temp_c', 'station_precip_mm', 'calc_ndvi_avg',
'calc_ndvi_range', 'calc_temp_range_c', 'calc_total_avg_air_temp',
'calc_avg_precip'],
dtype='object')
Como los algoritmos a usar están basado en el cálculo de distancias, primero se tendrán que normalizar los datos, de forma que todas las distancias sean relativas en una misma escala (entre 0 y 1).
Normalización: MinMaxScaler¶
#http://scikit-learn.org/stable/modules/preprocessing.html
from sklearn.preprocessing import MinMaxScaler
min_max_scaler = MinMaxScaler()
features_norm = min_max_scaler.fit_transform(train_transposed)
Matriz de similitud¶
Como primer acercamiento al agrupamiento entre características usando métodos de clasificación, se usará la matriz de similitud. Esta matriz utiliza la distancia entre instancias (que en este nuevo dataset traspuesto, las instancias son las características) para agrupar según su cercanía.
Además de calcular las distancias,
#http://docs.scipy.org/doc/scipy/reference/cluster.html
from scipy import cluster
from sklearn import metrics
dist = metrics.DistanceMetric.get_metric('euclidean')
matdist = dist.pairwise(features_norm)
Ahora tengo una matriz cuadrada de la distancia entre las caraceterísticas. Para plotearla de forma que estas distancias queden normalizadas entre 0 y 1, se usará de nuevo el MinMaxScaler. Sin embargo, en esta ocasión, se realizará un cambio importante.
Esta función realiza un escalado por columna en un dataframe o numpy array. Sin embargo, para nuestra matriz de distancias, nos interesa que el escalado se haga de forma uniforme para toda la matriz sin tener en cuenta a la columna que pertenece el dato, ya que estamos realizando comparaciones de distancias generales.
Para ello, se "aplanará" (flatten()) la matriz cuadrada, se redimensionará de forma que todos los datos estén en una sola columna, y por último, se redimensionará de nuevo a la forma original.
matdist_scaled = min_max_scaler.fit_transform(matdist.flatten().reshape(-1, 1))
matdist_scaled = matdist_scaled.reshape(matdist.shape[0], matdist.shape[1])
Ahora sí, se plotea el heatmap
#Se usa una función definida anteriormente para mantener formato en heatmaps
plot_heatmap(matdist_scaled, names, max = 1, min = 0)#Se usa una función definida anteriormente para mantener formato en heatmaps
Como se puede observar, existen muchas variables cuya distancia es cero sin ser el cruce de dicha característica consigo misma. Esto es debido a que, aunque su valor real no es cero, es muy cercano a cero. Al escalarlo entre 0 y 1, aún se acerca más al cero, y solo se muestran dos decimales.
Por lo tanto, este heatmap es una muy buena aproximación para visualizar los grupos de caraceterísticas que están muy cerca los unos de los otros en cuanto a distancia de valores se refiere.
Por ejemplo, todo el grupo de reanalysis_xxx_temp_kestán muy cerca entre sí, pero muy alejados del resto de características. Esto es un buen indicador de la presencia de un cluster bien diferenciado.
Jerarquía mediante Linkage¶
Otra manera de visualizar posibles formas de agrupar las características, es usando el algoritmo jerárquico.
La forma en la que funciona este algoritmo realizando sus cálculos será basándose en la matriz de similitud vista anteriormente. Sin embago, esta jerarquización y unión de instancias (características en este caso) dependerá del método que se use en el algoritmo.
En primer lugar, se usará el método de cálculo de distancias single, que toma la mínima distancia de las instancias que forman un cluster al compararlo con otro cluster para formar nuevos clusteres mayores.
Además, se visualizará mediante un dendograma.
from scipy import cluster
# http://docs.scipy.org/doc/scipy/reference/generated/scipy.cluster.hierarchy.linkage.html#scipy.cluster.hierarchy.linkage
clusters = cluster.hierarchy.linkage(matdist, method = 'single')
# http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.cluster.hierarchy.dendrogram.html
k = 2
fig, axs = plt.subplots(k, k, figsize = (10, 10))
#Utilizo el generador creado anteriormente para subplots cuadrados
incr_i_j = increment_index(k)
for cut in (5, 10, 15, 20):
i, j = next(incr_i_j)
cluster.hierarchy.dendrogram(clusters, color_threshold = cut, leaf_font_size = 8,
labels = names , leaf_rotation=90, ax = axs[i, j])
axs[i, j].set_title(f'Cut: {cut}', fontsize = 8)
plt.show()
<ipython-input-66-8619e0d2c8ac>:2: ClusterWarning: The symmetric non-negative hollow observation matrix looks suspiciously like an uncondensed distance matrix clusters = cluster.hierarchy.linkage(matdist, method = 'single')
Se realizan varios cortes sobre los umbrales (5, 10, 15, 20), obteniendo diferentes agrupaciones. Estos cortes representan el valor de distancias a partir del cual se dejan de formar nuevos clusteres.
De forma generalizada, cuanto mayor es la distancia entre lineas horizontales, mayor es la distancia entre los nuevos clusteres y, por lo tanto, mejor sería la agrupación.
Sin embargo, no siempre tiene que ser así. En este ejemplo, si se toma el valor de 20, se observa que se forman solamente dos clusteres. Como era de esperar por las observaciones de la matriz de similitud, las variables reanalysis_xxx_air_temp_kes uno de esos grupos. De hecho, para cada iteración han formado su propio cluster.
Desde mi punto de vista, la mejor agrupación es aquella cuyo umbral está en 10. Este umbral define 4 grupos, de los cuales, 1 es una única característica. Se trata de reanalysis_relative_humidity_percenty se interpreta como un outlier.
Para visualizar esta agrupación en 2d, se va a hacer uso de una herramiento de reducción de dimensionalidad: Análisis de Componentes Principales (PCA).
PCA (Principal Component Analysis)¶
A continuación, y para poder representar en 2d los clusters de características, se utilizará el algoritmo PCA para reducir la dimensionalidad a 2.
from sklearn.decomposition import PCA
pca = PCA (n_components = 2)
X_pca = pca.fit_transform(features_norm)
print("Componentes lineales:\n", pca.components_)
print("\nRatio de variabilidad: ", pca.explained_variance_ratio_, "\n")
Componentes lineales: [[ 0.02673872 0.02674034 0.02662135 ... 0.02609486 0.02594988 0.02594301] [-0.01404983 -0.01882628 -0.00384252 ... -0.00368385 0.00231955 -0.00321332]] Ratio de variabilidad: [0.97857145 0.0132881 ]
Como se puede observar, PCA ha hecho un gran trabajo, pues con tan solo dos componentes ha sido capaz de recoger casi el 98% de la varianza de las características originales.
Se realiza un corte por el umbral definido
cut = 10
labels = cluster.hierarchy.fcluster(clusters, cut , criterion = 'distance')
labels
array([2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 3, 4, 3, 2, 2, 2, 2, 2, 2, 3, 2, 2,
2, 2, 3], dtype=int32)
Se plotea en 2D los puntos proyectados de la PCA.
plt.scatter(X_pca[:,0], X_pca[:,1], c=labels, s=50, cmap="Set2")
for i in range(len(X_pca)):
plt.text(X_pca[i][0], X_pca[i][1], names[i], fontsize = 8)
plt.show()
Gracias a las proyecciones PCA, se puede visualizar claramente los distintos clusteres que se han realizado en base al algoritmo jerárquico. De hecho, tienen bastante sentido visual (esto podría no ser así aun teniendo sentido matemático).
Se comprueba que reanalysis_relative_humidity_percentes una característica outlier. Esto, en este contexto, hace que dicha caraceterística sea muy valiosa, ya que en caso de necesitar reducir dimensiones eliminando variables, esta característica sea dificilmente reemplazable.
De las demas, se podrían realizar técnicas de reducción de la dimensionalidad por clusteres para quedarnos con las más representativas (por ejemplo, un PCA por cluster).
Esto se haría en caso de que computacionalmente el modelo sea demasiado pesado.
Como no se visualizan muy bien las características agrupadas en la esquina inferior izquierda, se va a representar una gráfica haciendo zoom en esas caraceterísticas.
xmax = 0
xmin = -10
ymax = -0.5
ymin = -1
for i in range(len(X_pca)):
if (X_pca[i][0] > xmin and
X_pca[i][1] < xmax and
X_pca[i][1] > ymin and
X_pca[i][1] < ymax):
plt.scatter(X_pca[i][0], X_pca[i][1], c=labels[i],s=50, cmap="Set2")
plt.text(X_pca[i][0], X_pca[i][1], names[i], fontsize = 8)
plt.show()
Con esto se da por finalizado la fase de ingeniería de características.
K-MEANS¶
Como primer algoritmo de clusterización para los datos se va a usar K-means
from sklearn.cluster import KMeans
De nuevo, se han de normalizar los datos, pues antes se normalizaó el dataset traspuesto.
train
| city | weekofyear | ndvi_ne | ndvi_nw | ndvi_se | ndvi_sw | reanalysis_air_temp_k | reanalysis_avg_temp_k | reanalysis_max_air_temp_k | reanalysis_min_air_temp_k | ... | station_avg_temp_c | station_diur_temp_rng_c | station_max_temp_c | station_min_temp_c | station_precip_mm | calc_ndvi_avg | calc_ndvi_range | calc_temp_range_c | calc_total_avg_air_temp | calc_avg_precip | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 18 | 0.122600 | 0.103725 | 0.198483 | 0.177617 | 297.572857 | 297.742857 | 299.8 | 295.9 | ... | 25.442857 | 6.900000 | 29.4 | 20.0 | 16.0 | 0.150606 | 0.094758 | 9.4 | 24.932857 | 24.000 |
| 1 | 1 | 19 | 0.169900 | 0.142175 | 0.162357 | 0.155486 | 298.211429 | 298.442857 | 300.9 | 296.4 | ... | 26.714286 | 6.371429 | 31.7 | 22.2 | 8.6 | 0.157479 | 0.027725 | 9.5 | 25.887857 | 13.270 |
| 2 | 1 | 20 | 0.032250 | 0.172967 | 0.157200 | 0.170843 | 298.781429 | 298.878571 | 300.5 | 297.3 | ... | 26.714286 | 6.485714 | 32.2 | 22.8 | 41.4 | 0.133315 | 0.140717 | 9.4 | 26.172857 | 33.750 |
| 3 | 1 | 21 | 0.128633 | 0.245067 | 0.227557 | 0.235886 | 298.987143 | 299.228571 | 301.4 | 297.0 | ... | 27.471429 | 6.771429 | 33.3 | 23.3 | 4.0 | 0.209286 | 0.116433 | 10.0 | 26.654286 | 8.950 |
| 4 | 1 | 22 | 0.196200 | 0.262200 | 0.251200 | 0.247340 | 299.518571 | 299.664286 | 301.9 | 297.5 | ... | 28.942857 | 9.371429 | 35.0 | 23.9 | 5.8 | 0.239235 | 0.066000 | 11.1 | 27.655714 | 9.000 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 1451 | 0 | 21 | 0.342750 | 0.318900 | 0.256343 | 0.292514 | 299.334286 | 300.771429 | 309.7 | 294.5 | ... | 28.633333 | 11.933333 | 35.4 | 22.4 | 27.0 | 0.302627 | 0.086407 | 13.0 | 27.408810 | 36.000 |
| 1452 | 0 | 22 | 0.160157 | 0.160371 | 0.136043 | 0.225657 | 298.330000 | 299.392857 | 308.5 | 291.9 | ... | 27.433333 | 10.500000 | 34.7 | 21.7 | 36.6 | 0.170557 | 0.089614 | 13.0 | 26.306667 | 121.850 |
| 1453 | 0 | 23 | 0.247057 | 0.146057 | 0.250357 | 0.233714 | 296.598571 | 297.592857 | 305.5 | 292.4 | ... | 24.400000 | 6.900000 | 32.2 | 19.2 | 7.4 | 0.219296 | 0.104300 | 13.0 | 23.924286 | 29.000 |
| 1454 | 0 | 24 | 0.333914 | 0.245771 | 0.278886 | 0.325486 | 296.345714 | 297.521429 | 306.1 | 291.9 | ... | 25.433333 | 8.733333 | 31.2 | 21.0 | 16.0 | 0.296014 | 0.088143 | 10.2 | 24.314524 | 39.165 |
| 1455 | 0 | 25 | 0.298186 | 0.232971 | 0.274214 | 0.315757 | 298.097143 | 299.835714 | 307.8 | 292.3 | ... | 27.475000 | 9.900000 | 33.7 | 22.2 | 20.4 | 0.280282 | 0.082786 | 11.5 | 26.211071 | 28.650 |
1446 rows × 25 columns
train_norm = min_max_scaler.fit_transform(train)
Como, a priori, se desconoce el número de clusters de los que se debería de componer el dataset, se realizará un primer análisis variando la cantidad de clusteres y comparando dos indicadores de cuan bueno es dicho cluster:
Distorsión: se trata de la suma de las distancias cuadradas medidas desde los centroides de cada cluster a todos los puntos de su cluster. Se busca la menor distorsión posible sin excederse en el número de clusters. Como caso ideal de resultado de distorsión igual a cero, sería tantos clusteres como puntos haya. Sin embargo, esta clusterización no tienen ningún sentido, ya que no se habría agrupado absolutamente nada. Para alcanzar este equilibrio se hará uso del método no cientítico del codo, pero que da muy buenas aproximaciones experimentales.
Silueta: este indicador da una medida promedio de dos conceptos fundamentales en la medición de cuán bueno es un cluster formado.
- Cohesión: cómo de cerca esta cada punto de un cluster con respecto al resto de puntos de ese mismo cluster.
- Separación: cómo de lejos está cada punto de un cluster con respecto a cada punto de los ostros clústeres.
El valor de la Silueta puede variar desde -1 (punto categorizado erróneamente) hasta 1 (punto perfectamente agrupado). Los valores alrededor de 0 son puntos justo en la frontera entre clusteres. Como la silueta da valores promedios de este cálculo entre todos los puntos, valores a 0 se interpretan como superposición de clusteres (es decir, mala definición de los clusteres).
distortions = []
silhouettes = []
for i in range(2, 20):
km = KMeans(i, init='k-means++', n_init=10, random_state=42)
clustering = km.fit_predict(train_norm)
distortions.append(km.inertia_)
silhouettes.append(metrics.silhouette_score(train_norm, clustering))
plt.plot(range(2,20), distortions, marker='o')
plt.xticks(range(2, 20))
plt.xlabel('K')
plt.ylabel('Distortion')
plt.grid(True)
plt.show()
Como era de esperar, la distorsión disminuye a medida que se crean más clusteres. Como se ha comentado anteriormente, es difícil determinar un buen punto donde hay un balance entre número de clusters y baja distorsión. El método del codo nos dice que justo en el punto donde hay un cambio de pendiente, se encuentra este balance.
Ojo, esto no es un cálculo 100% efectivo, solo orientativo. El problema es que al haber tantas dimensiones de caraceterísticas no es posible aplicar el sentido común para escoger este número de clústeres.
Por lo tanto, se usará este método auxiliar.
!pip install kneed
import kneed
Collecting kneed Downloading kneed-0.8.5-py3-none-any.whl.metadata (5.5 kB) Requirement already satisfied: numpy>=1.14.2 in /usr/local/lib/python3.10/dist-packages (from kneed) (1.26.4) Requirement already satisfied: scipy>=1.0.0 in /usr/local/lib/python3.10/dist-packages (from kneed) (1.13.1) Downloading kneed-0.8.5-py3-none-any.whl (10 kB) Installing collected packages: kneed Successfully installed kneed-0.8.5
kneedle = kneed.KneeLocator(range(2, 20), distortions[:20], curve="convex",
direction="decreasing")
elbow_point = kneedle.elbow
print('Elbow: ', elbow_point)
kneedle.plot_knee()
Elbow: 6
Según el método del codo, la mejor agrupación sería con 6 clústeres.
Ahora, se va a plotear los resultados del indicador de silueta.
def plot_silhouettes(silhouettes, min, max, xLabel, step = 1):
plt.plot(range(min, max, step), silhouettes , marker='o')
plt.xticks(range(min, max, step))
plt.xlabel(xLabel)
plt.ylabel('Silhouette')
plt.tick_params(axis = 'x', labelsize = 8)
plt.grid(True)
plt.show()
print("Valores de silueta:")
print([round(silhouette, 2) for silhouette in silhouettes])
plot_silhouettes(silhouettes, 2, 20, 'K Clusters')
Valores de silueta: [0.47, 0.38, 0.31, 0.28, 0.28, 0.28, 0.22, 0.21, 0.21, 0.2, 0.2, 0.15, 0.16, 0.15, 0.15, 0.14, 0.14, 0.13]
Se obtiene una máxima silueta para un valor k = 2 clústeres.
Este resultado es contradictorio al del método del codo, que sugiere 6 clústeres.
Esto ocurre porque el método del codo no tiene en cuenta la separación de clústeres, tan solo la compactación en cada uno de ellos. Para el método del codo, cuantos más clústeres haya, más compactos estarán los datos intracluster, hasta que llega un momento (punto de inflexión) en el que la mejora ya no es significativa.
La silueta, en cambio, además de tener en cuenta la compactación o cohesión, también tiene en cuenta la clara separación entre clústeres.
Por lo tanto, se considera un mejor medidor de la calidad del modelo de clusterización.
Por otra parte, el mejor resultado obtenido con 2 clústeres es de 0.47. Este valor, cercano a 0.5, sugiere una moderada/buena agrupación. Quizás haya un porcentaje importante de puntos por los bordes y esto hace que la silueta no sea mayor, pero no es un valor cercano a 0, que sí sería un mal resultado.
No obstante, se va a visualizar cómo quedaría la agrupación con 2 y 6 clústeres respectivamente en 2d usando PCA.
X_pca = pca.fit_transform(train_norm)
Se define una función para usarse recurrentemente (por este y más modelos) para plotear los puntos de train_norm en 2D
def plot_train_norm_2D(clustering, algorithm_type, size = 10):
plt.scatter(X_pca[:, 0], X_pca[:, 1], s=size, c = clustering, cmap="Set3",
zorder = 0)
plt.title(f"Clustering con {algorithm_type} en 2D")
plt.xlabel("Componente Principal 1")
plt.ylabel("Componente Principal 2")
plt.grid(True)
plt.show()
#Referencia: apuntes de clase
Se calcula la clusterización de kmeans para 2 clústeres. Además, se plotean antes de llamar a la función antes definida.
km = KMeans(2, init = 'k-means++', n_init = 10, random_state = 42)
clustering = km.fit_predict(train_norm)
centroids_2d = pca.fit_transform(km.cluster_centers_)
plt.scatter(centroids_2d[:, 0], centroids_2d[:, 1], marker='x', s=100, c='red')
plot_train_norm_2D(clustering, 'K-Means')
km = KMeans(elbow_point, init = 'k-means++', n_init = 10, random_state = 42)
clustering = km.fit_predict(train_norm)
centroids_2d = pca.fit_transform(km.cluster_centers_)
plt.scatter(centroids_2d[:, 0], centroids_2d[:, 1], marker='x', s=100, c='red')
plot_train_norm_2D(clustering, 'K-Means')
Usando la visualización, parece que la subdivisión en 2 clústeres tiene mucho sentido. Sin embargo, no hay que olvidar que lo que se ve es una proyección de muchas dimensiones a 2 y esto hace que se dificulte enormemente el análisis visual.
DBSCAN¶
DBSCAN es un algoritmo de agrupamiento basado en la densidad de puntos. Otra gran diferencia con respecto a k-means, es que no es necesario parametrizar, a priori, el número de clústeres que se desean formar.
En primer lugar, se importa la librería.
from sklearn.cluster import DBSCAN
Como un primer acercamiento, se va a hacer igual que se hizo con kmeans: varíando uno de sus parámetros clave, vamos a observar como varían un par de parámetros de control de calidad de la dispersión.
En este caso, se usará igualmente la silueta pero, además, el porcentaje de puntos de ruido. Estos puntos de ruido son puntos no clasificados en ningún cluster, es decir, puntos que se han quedado fuera de épsilon para todos los épsilon posibles formados.
De manera general, se suelen tomar como puntos mínimos, al menos, el número de dimensiones + 1. Como esto es solo orientativo, se va a iterar desde 2 hasta 40 puntos mínimos de agrupamiento.
silhouettes = []
porc_noise_points = []
for i in range(2, 40):
dbscan = DBSCAN(eps=0.5, min_samples = i)
clustering = dbscan.fit_predict(train_norm)
silhouettes.append(metrics.silhouette_score(train_norm, clustering))
noise_points_index = np.where(clustering == -1)[0]
noise_points = train_norm[noise_points_index]
porc_n_p = round(noise_points.shape[0] / train_norm.shape[0], 2)
porc_noise_points.append(porc_n_p)
plot_silhouettes(silhouettes, 2, 40, 'min_samples')
Valores de silueta: [0.11, 0.33, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.39, 0.39, 0.39, 0.38, 0.38, 0.38, 0.38, 0.37, 0.37, 0.37, 0.37, 0.37, 0.37, 0.36, 0.36, 0.36, 0.36, 0.36, 0.35, 0.35, 0.35, 0.34, 0.34, 0.34, 0.33, 0.33, 0.33, 0.33]
Se grafican, también, los puntos de ruido en función del min_samples
plt.plot(range(2,40), porc_noise_points , marker='o')
plt.xticks(range(2, 40))
plt.xlabel('min_samples')
plt.ylabel('% Puntos de Ruido')
plt.tick_params(axis = 'x', labelsize = 5)
plt.grid(True)
plt.show()
Los resultados obtenidos no malos, pero ligeramente peores que con K-Means. El máximo puntaje de silueta está un poco por encima de 0.40 para 4 clusteres.
Sin embargo, estas gráficas muestran la evolución del modelo variando tan solo un hiperparámetro, el MinPts (número mínimo de puntos de muestra que requiere un punto para ser considerado como centro para formar un cluster). Al ser un modelo que requiere de dos parámetros principales, épsilon y MinPts (aunque tiene más, estos son los más influyentes), se va a realizar una exploración de la combinatoria de ambas variables.
Para esta tarea se puede pensar en el uso de funciones especializadas para ello como GridSearchCV. Esta función combina todas las posibilidades que se configuren y se quedará con la mejor según el tipo de evaluación que se le de.
Sin embargo, esta función no integra la métrica de la silueta, ya que todas las métricas que incorpora son de aprendizaje supervisado, es decir, esperan tener etiquetas (como métricas de homogeneidad y completitud).
Por lo tanto, se creará una función específica para esta labor con la métrica de la silueta.
Primero, Se definen aquellos hiperparámetros que se quieren explorar, así como sus rangos de exploración (igual que se haría con RandomGridSearchCV).
param_grid = {
'eps': np.arange(0.2, 1, 0.1),
'min_samples': range(2, 40)
}
A continuación, se crea el algoritmo que itera entre todos los valores de los dos hiperparámetros y sus combinaciones. Si el resultado es mejor que el anterior, sobreescribe el mejor resultado y al final
best_sil_score = 0
best_params = {key: 0 for key in param_grid.keys()}
for value1 in param_grid['eps']:
for value2 in param_grid['min_samples']:
dbscan = DBSCAN(eps=value1, min_samples = value2)
clustering = dbscan.fit_predict(train_norm)
n_labels = len(np.unique(clustering))
if n_labels == 1: #Esto se hace para evitar un error de silhouete_score
continue
sil_score = metrics.silhouette_score(train_norm, clustering)
if sil_score > best_sil_score:
best_sil_score = sil_score
best_params['eps'] = value1
best_params['min_samples'] = value2
noise_points_index = np.where(clustering == -1)[0]
noise_points = train_norm[noise_points_index]
porc_n_p = round(noise_points.shape[0] / train_norm.shape[0], 2)
print("Mejores hiperparámetros:", best_params)
print("Mejor puntuación de la silueta:", best_sil_score)
print("Ratio de puntos de ruido:", porc_n_p)
Mejores hiperparámetros: {'eps': 0.9000000000000001, 'min_samples': 2}
Mejor puntuación de la silueta: 0.4671721515440993
Ratio de puntos de ruido: 0.0
Veamos el número de clusters formados con estos hiperparámetros
dbscan = DBSCAN(eps=0.9, min_samples = 2)
clustering = dbscan.fit_predict(train_norm)
np.unique(clustering)
array([0, 1])
Se han formado dos clusters:
- Grupo 0
- Grupo 1
No hay grupo -1, es decir, no hay puntos de ruido.
Vamos a representar gráficamente estos grupos en 2d con PCA.
plot_train_norm_2D(clustering, 'DBSCAN')
Tanto la silueta como la inspección visual indican que esta clusterización coincide con la mejor de K-means. Esto ocurre porque, aunque son dos algoritmos distintos, tienen un punto en común:
- Ambos algoritmos detectan clústeres globulares (un círculo, esfera, hiperesfera...).
Esto no quiere decir necesariamente que el mejor tipo de forma de cluster para este dataset sea en forma de hiperesfera, pero al menos se sabe que esta forma ha aglomerado de forma similar en dos tipo de algoritmos matemáticamente diferentes.
Además, el resultado es aceptable.
EXPECTATION - MAXIMIZATION: Gaussian Mixture Models¶
Para probar un enfoque distinto a los anteriores se harán pruebas con los modelos mixtos gaussianos.
Este tipo de algoritmo, al contrario que los anteriores, utiliza un método generativo en lugar de discriminativo. Es decir, se calcula el modelo entrando con los datos del dataset y, una vez creado, es capaz de calcular la probabilidad de cada punto de pertenecer a un cluster e incluso crear su propia nube de puntos en función de su función probabilística Gaussiana.
from sklearn.mixture import GaussianMixture
silhouettes = []
for i in range(2, 20):
gmm = GaussianMixture(n_components=i)
clustering = gmm.fit_predict(train_norm)
silhouettes.append(metrics.silhouette_score(train_norm, clustering))
plot_silhouettes(silhouettes, 2, 20, 'GMM')
Valores de silueta: [0.47, 0.35, 0.32, 0.23, 0.24, 0.16, 0.23, 0.11, 0.14, 0.16, 0.11, 0.11, 0.09, 0.12, 0.08, 0.09, 0.09, 0.08]
gmm = GaussianMixture(n_components=2)
clustering = gmm.fit_predict(train_norm)
plot_train_norm_2D(clustering, 'GMM')
metrics.silhouette_score(train_norm, clustering)
0.4671721515440993
Al ser GMM un modelo probabilístico también se puede encontrar la probabilidad de que un punto pertenezca al cluster asignado.
Esto lo podemos hacer usando la función predict_proba, que devuelve una matriz de tamaño [n_samples, n_clusters] que mide la probabilidad de cada punto de pertenecer a cada cluster.
probs = gmm.predict_proba(train_norm)
probs[:30]
array([[0., 1.],
[0., 1.],
[0., 1.],
[0., 1.],
[0., 1.],
[0., 1.],
[0., 1.],
[0., 1.],
[0., 1.],
[0., 1.],
[0., 1.],
[0., 1.],
[0., 1.],
[0., 1.],
[0., 1.],
[0., 1.],
[0., 1.],
[0., 1.],
[0., 1.],
[0., 1.],
[0., 1.],
[0., 1.],
[0., 1.],
[0., 1.],
[0., 1.],
[0., 1.],
[0., 1.],
[0., 1.],
[0., 1.],
[0., 1.]])
Visualizando los primeros 30 puntos, todos pertenecen a su clase con una probabilidad del 100%. ¿Se cumplirá esto para todos los puntos del dataset? Veamos los valores únicos de probabilidades.
np.unique(probs)
array([0., 1.])
Efectivamente, de acuerdo a este tipo de clasificación la probabilidad es totalmente determinante.
Parece cada vez más claro que el número de clústeres óptimo para este dataset es 2.
Son 3 los modelos que han agrupado exactamente igual, teniendo los 3 distintas bases matemáticas y, para este en concreto, una forma geométrica de clusterizar diferente a las dos anteriores: en forma de elipse
Visualizemos, además, las elipses que conforman la definición de la gaussiana y sus curvas de nivel.
from matplotlib.patches import Ellipse
def draw_ellipse(position, covariance, ax=None, **kwargs):
"""Draw an ellipse with a given position and covariance"""
ax = ax or plt.gca()
# Convert covariance to principal axes
if covariance.shape == (2, 2):
U, s, Vt = np.linalg.svd(covariance)
angle = np.degrees(np.arctan2(U[1, 0], U[0, 0]))
width, height = 2 * np.sqrt(s)
else:
angle = 0
width, height = 2 * np.sqrt(covariance)
# Draw the Ellipse
for nsig in range(1, 4):
ax.add_patch(Ellipse(position, nsig * width, nsig * height,
angle, **kwargs))
def plot_gmm(gmm, X, label=True, ax=None):
ax = ax or plt.gca()
labels = gmm.fit(X).predict(X)
if label:
ax.scatter(X[:, 0], X[:, 1], c=labels, s=10, cmap='viridis', zorder=2)
else:
ax.scatter(X[:, 0], X[:, 1], s=10, zorder=2)
ax.axis('equal')
w_factor = 0.2 / gmm.weights_.max()
for pos, covar, w in zip(gmm.means_, gmm.covariances_, gmm.weights_):
draw_ellipse(pos, covar, alpha=w * w_factor)
ax.scatter(gmm.means_[:, 0], gmm.means_[:, 1], marker='*', c='r', s=300, zorder=3)
#Referencia: apuntes aportados en clase
plot_gmm(gmm, X_pca)
plt.grid(True)
plt.show()
<ipython-input-99-7ca1e15ecad3>:18: MatplotlibDeprecationWarning: Passing the angle parameter of __init__() positionally is deprecated since Matplotlib 3.6; the parameter will become keyword-only two minor releases later. ax.add_patch(Ellipse(position, nsig * width, nsig * height,
Para modelos probabilísticos, se pueden aplicar otras métricas para comprobar cuán bueno es dicho modelo. Una de estas métricas es el (AIC)](https://en.wikipedia.org/wiki/Akaike_information_criterion). Con esta métrica se hará una comparación de la cantidad de pérdida de información en función del modelo.
Por lo tanto, se realizará esta prueba variando el número de componentes (clústeres).
n_components = np.arange(2, 20)
models = [GaussianMixture(n, covariance_type='full', random_state=42)
for n in n_components]
aics = [model.fit(train_norm).aic(train_norm) for model in models]
plt.plot(n_components, aics, marker = 'o')
plt.grid(True)
#Referencia: apuntes de clase
Se extrae un número de componentes óptimos para minimizar el AIC de 8 componentes. Este es un número de clústeres que nada tiene que ver con lo calculado anteriormente usando la silueta. Veámoslo.
gmm = GaussianMixture(n_components=8)
clustering = gmm.fit_predict(train_norm)
plot_train_norm_2D(clustering, 'GMM')
metrics.silhouette_score(train_norm, clustering)
0.20193210651635332
Ya no es tan intuitivo al visualizarlo en 2D. Pero qué es más fiable, ¿El AIC o la Silueta?
Realmente, ninguno es mejor que el otro. Son dos puntos de vista de comprobar cuán bueno es un modelo.
La silueta es una medida de clusterización, y busca clústeres separados y bien definidos.
AIC se centra en cuán bueno es GMM como estimador de densidad, es decir, prioriza un buen ajuste estadístico. Esto favorece un ajuste debido a cualquier mínima variación en los datos, lo que lleva a una mayor fragmentación de los clústeres aunque no queden separados y visualmente bien divididos. Obviamente, esto va en detrimento de la Silueta.
Si se busca tener una clara definición de grupos, la Silueta sería la mejor opción como medida y 2 grupos sería lo ideal.
Sin embargo, si se busca el realizar un etiquetado de datos en base a un modelo que capture las variaciones internas de los datos para que, con estas etiquetas, se pueda realizar una predicción, el AIC sería la medida a escoger y, por lo tanto, 8 clústeres sería lo ideal.
Otra posibilidad, es encontrar un punto adecuado de encuentro de ambas métricas.
fig, ax1 = plt.subplots()
ax1.plot(range(2, 20), silhouettes , marker='o',
color = 'salmon', label = 'Silhouette')
ax1.set_xticks(range(2, 20))
ax1.set_ylabel('Silhouette')
ax2 = ax1.twinx()
ax2.plot(n_components, aics, marker = 'o',
color = 'skyblue', label = 'AIC')
ax2.set_ylabel('AIC')
ax1.legend(loc = 'upper right')
ax2.legend(loc = 'upper left')
plt.axvline(x = 3, color = 'black', linestyle = '--')
plt.grid(True)
plt.tight_layout()
plt.show()
En esta gráfica se ve claramente como con 2 clústeres la Silueta es aceptable pero el AIC es relativamente malo, y con 8 clústeres precisamente la silueta cae por debajo de 0.20.
Sin embargo, para 3 clústeres, el AIC cae drásticamente con respecto a 2 (aunque no alcanza sus valores más bajos) y al AIC se mantiene cercano a 0.35, que sigue siendo un valor no bueno, pero aceptable.
Esto puede ser un punto medio de encuentro entre un modelo estadísticamente ajustado y unas agrupaciones que, en este punto de la silueta son mediocres, pero son algo mejores que en el mejor punto del AIC.
Sin embargo, es necesario considerar cual sería el objetivo de la clusterización para optar, en este caso, por el número de clústeres.
AFFINITY PROPAGATION¶
AP es un algoritmo de agrupamiento que no necesita la especificación previa del número de clústeres, al igual que como pasaba con DBSCAN. Sin embargo, a diferencia de este modelo, AP selecciona automáticamente "ejemplares" de la matriz de similitud para ser centros de clústeres.
El parámetro que controla la cantidad de estos puntos ejemplares es preference. Se suele escoger un número negativo ya que, a menor valor de este hiperparámetro, menos puntos serán escogidos como posibles centros, y menor será la cantidad de clusteres (aunque por aleatoriedad esto no se cumple a rajatabla).
Como se ha hecho en los otros algoritmos, se va a realizar un estudio de la silueta en función de su parámetro principal, en este caso, preference.
from sklearn.cluster import AffinityPropagation
silhouettes = []
for i in range(-110, 0, 10):
af = AffinityPropagation(preference = i, random_state=42,
max_iter = 1000, convergence_iter = 30)
clustering = af.fit_predict(train_norm)
if len(np.unique(clustering)) == train_norm.shape[0]:
silhouettes.append(0)
else:
silhouettes.append(metrics.silhouette_score(train_norm, clustering))
/usr/local/lib/python3.10/dist-packages/sklearn/cluster/_affinity_propagation.py:142: ConvergenceWarning: Affinity propagation did not converge, this model may return degenerate cluster centers and labels. warnings.warn(
plot_silhouettes(silhouettes, -110, 0, 'Affinity', 10)
Valores de silueta: [0.38, 0, 0.31, 0.38, 0.31, 0.29, 0.31, 0.31, 0.28, 0.26, 0.16]
Se observa que el mejor valor de preference es de -80. Por lo tanto, es el que se va a usar para generar el modelo.
Se observa que el valor de la silueta es menor que en los modelos anteriores, por lo que se espera un resultado diferente en la clusterización.
af = AffinityPropagation(preference = -80, random_state=42,
max_iter = 1000, convergence_iter = 30)
clustering = af.fit_predict(train_norm)
Se comprueba la cantidad de clústeres:
np.unique(clustering)
array([0, 1, 2])
plot_train_norm_2D(clustering, 'AP')
metrics.silhouette_score(train_norm, clustering)
0.37815382611842574
Este modelo devuelve una asociación de 3 clústeres para el dataset. Es lógico que se obtenga este resultado diferente donde, la silueta, tiene un valor más bajo.
Esto no es necesariamente un punto en contra.
Este algoritmo en lugar de trabajar con las distancias euclideas directamente, lo hace con la matriz de similitud (igual que los algoritmos jerárquicos). Por lo tanto, puede ser que la asociación tenga más en cuenta similitudes de puntos que su distancia directa medida punto a punto.
Los resultados que se han obtenido en los cuatro modelos, de forma general, ha sido mediante la métrica de la Silueta.
Estos resultados han sido similares en todos ellos: una silueta con un valor de 0.47 y un número óptimo de 2 clusters.
Este valor de silueta no es muy bueno, pero es aceptable. Además, en las visualizaciones en 2D usando PCA sobre los datos normalizados, se demuestra que 2 clústeres parece lo más razonable.
Sin embargo, según el modelo probabilístico de Expectation-Maximization Gaussian Mixture Models y la prueba del AIC, se demuestra que con 8 clústeres se minimiza la pérdida de información, obteniendo un modelo más sensible a variaciones de los datos.
El tener más clústeres estándo creando un modelo más ajustado probabilísticamente pero con fronteras menos claras, puede tener la ventaja que, al etiquetar estos valores, puede que se cree una nueva característica que sirva como una dimensión más para posibles futuras predicciones que, a diferencia de un etiquetado de dos clústeres, tenga un nivel de profundización en la explicación del modelo mucho mayor. Sin embargo, pierde claridad de clusterización al no haber grupos claramente separados y añade complejidad a un futuro modelo donde dichas etiquetas se convertirian en múltiples columnas, ya que se tendría que realizar un OneHotEncoding (al no ser una variable ordinal).
Por lo tanto, se va a escoger un punto intermedio también hayado mientras se exploraba el modelo Gaussian Mixture Models, donde se tienen 3 clústeres. Esta elección viene, además, reforzada por el resultado obtenido usando el modelo de Affinity Propagation (aunque la agrupación en sí es distinta).
Con 3 clústeres, la métrica de la silueta muestra una caída de un punto con respecto a una clusterización de 2 grupos, pero bastante más alto que con 8 clústeres.
La razón por la que se sacrifica claridad de separación entre grupos en pos de la métrica del AIC, es porque se quiere aprovechar la capacidad del modelo de generar etiquetas que sirvan como buenas predictoras para un futuro trabajo de aprendizaje supervisado.
Además, el AIC es una métrica que penaliza la cantidad de características (clústeres en este caso), por lo que asegura que no se da overfitting ni underfitting.
Para este trabajo futuro, se podría reutilizar todo el análisis exploratorio de datos añadiendo una exploración concreta con la etiqueta dada, que será la cantidad de casos de contagio de Dengue.
Un posible trabajo futuro, sería la exploración con otros tipos de algoritmos de clasificación así como un análisis de entendimiento de por qué se generan los clústeres que se han generado, etiquetando el dataset según estos clústeres, e ir explorando de nuevo todas las variables con respecto a este etiquetado como si de una variable objetivo se tratase.
A continuación y, para terminar, se muestra cómo quedarían los datos divididos en 3 clústeres con el modelo GMM y se exportará el dataset ya etiquetado con su cluster para que sirva de entrada en el siguiente ejercicio: Entrenamiento Supervisado
gmm = GaussianMixture(n_components=3)
clustering = gmm.fit_predict(train_norm)
plot_train_norm_2D(clustering, 'GMM')
metrics.silhouette_score(train_norm, clustering)
0.3522048606563886
Se etiqueta el dataset y se guarda con sus índices originales (para poder cruzar datos a posteriori si es necesario, ya que se han borrado columnas y filas)
train_norm = pd.DataFrame(train_norm, columns = train.columns)
train_norm['cluster'] = clustering
train_norm.to_csv('train_norm_clustered.csv', index = True)
train_norm.head()
| city | weekofyear | ndvi_ne | ndvi_nw | ndvi_se | ndvi_sw | reanalysis_air_temp_k | reanalysis_avg_temp_k | reanalysis_max_air_temp_k | reanalysis_min_air_temp_k | ... | station_diur_temp_rng_c | station_max_temp_c | station_min_temp_c | station_precip_mm | calc_ndvi_avg | calc_ndvi_range | calc_temp_range_c | calc_total_avg_air_temp | calc_avg_precip | cluster | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1.0 | 0.333333 | 0.578226 | 0.614835 | 0.386418 | 0.395544 | 0.388291 | 0.354667 | 0.123457 | 0.692308 | ... | 0.210393 | 0.174194 | 0.486239 | 0.052305 | 0.407631 | 0.141899 | 0.261745 | 0.358874 | 0.099174 | 2 |
| 1 | 1.0 | 0.352941 | 0.629943 | 0.657063 | 0.321190 | 0.359233 | 0.472710 | 0.441778 | 0.191358 | 0.730769 | ... | 0.163498 | 0.322581 | 0.688073 | 0.028114 | 0.419152 | 0.032783 | 0.268456 | 0.504184 | 0.054835 | 2 |
| 2 | 1.0 | 0.372549 | 0.479441 | 0.690881 | 0.311879 | 0.384430 | 0.548064 | 0.496000 | 0.166667 | 0.800000 | ... | 0.173638 | 0.354839 | 0.743119 | 0.135338 | 0.378645 | 0.216710 | 0.261745 | 0.547549 | 0.139463 | 1 |
| 3 | 1.0 | 0.392157 | 0.584823 | 0.770066 | 0.438912 | 0.491150 | 0.575260 | 0.539556 | 0.222222 | 0.776923 | ... | 0.198986 | 0.425806 | 0.788991 | 0.013076 | 0.505996 | 0.177182 | 0.302013 | 0.620802 | 0.036983 | 1 |
| 4 | 1.0 | 0.411765 | 0.658698 | 0.788882 | 0.481601 | 0.509943 | 0.645515 | 0.593778 | 0.253086 | 0.815385 | ... | 0.429658 | 0.535484 | 0.844037 | 0.018960 | 0.556200 | 0.095086 | 0.375839 | 0.773177 | 0.037190 | 1 |
5 rows × 26 columns
Para descargar el archivo en local
files.download('train_norm_clustered.csv')